#pragma once
extern "C" {
#include <f2c.h>
#include <clapack.h>
}
#include "IComplexMatrix.h"

namespace LatoolNet {
	using namespace System;

	ref class ComplexGeneralMatrix : public IComplexMatrix {
	internal:
		doublecomplex * dat;
	private:
		long * ipiv;
		int m_rownum;
		int m_colnum;
		bool m_isFactorized;

		ComplexGeneralMatrix(ComplexGeneralMatrix ^ orig) {
			m_rownum = orig->m_rownum;
			m_colnum = orig->m_colnum;
			dat = new doublecomplex [orig->m_rownum * orig->m_colnum];

			for (int i = 0; i < m_rownum * m_colnum; i++) {
				dat[i] = orig->dat[i];
			}
			ipiv = new long[m_rownum];

			for (int i = 0; i < m_rownum; i++) {
				ipiv[i] = orig->ipiv[i];
			}
			m_isFactorized = orig->m_isFactorized;

		};
	public:
		ComplexGeneralMatrix(int rownum, int colnum) {
			m_rownum = rownum;
			m_colnum = colnum;
			dat = new doublecomplex[rownum * colnum];

			for (int i = 0; i < rownum * colnum; i++) {
				dat[i].r = 0.0;
				dat[i].i = 0.0;
			}
			ipiv = new long[rownum];
			m_isFactorized = false;
		};

		~ComplexGeneralMatrix(){
			delete[] dat;
			delete[] ipiv;
		};

		virtual IComplexMatrix ^ Clone() {		
			return (IComplexMatrix^) gcnew ComplexGeneralMatrix(this);
		};

		virtual property int RowNum {
			int get() {
				return m_rownum;
			}
		};
		virtual property int ColNum {
			int get() {
				return m_colnum;
			}
		};
		virtual property Complex default[int, int] {
			Complex get(int i, int j) {
				return Complex(dat[j * m_rownum + i].r, dat[j * m_rownum + i].i);
			}
			void set(int i, int j, Complex value) {
				doublecomplex c;
				c.r = value.Real;
				c.i = value.Imag;
				dat[j * m_rownum + i] = c;
				m_isFactorized = false;
			}
		};
		virtual property MatrixType Type {
			MatrixType get() { return MatrixType::ComplexGeneral; }
		};
		virtual property bool IsFactorized {
			bool get() { return m_isFactorized; }
		};

		virtual void Invert() {

			long n = m_rownum;
			long lda = n;
			long lwork = n;

			doublecomplex * work = new doublecomplex[lwork];
			long info = 0;

			zgetri_(&n, dat, &lda, ipiv, work, &lwork, &info);

			delete[] work;

			m_isFactorized = false;
		};

		virtual void Factorize() {

			long m = m_rownum;
			long n = m_colnum;
			long lda = m;
			long info;

			zgetrf_(&m, &n, dat, &lda, ipiv, &info);

			m_isFactorized = true;
		};

		virtual void Solve(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			char trans = 'N';
			long n = m_rownum;
			long nrhs = gb->m_colnum;
			long lda = n;
			long ldb = n;
			long info = 0;

			zgetrs_(&trans, &n, &nrhs, dat, &lda, ipiv, gb->dat, &lda, &info);
		};

		virtual void SolveWithSimpleDriver(IMatrix ^ b) {

			ComplexGeneralMatrix ^ gb = (ComplexGeneralMatrix^) b;

			long n = m_rownum;
			long nrhs = gb->ColNum;
			long lda = n;
			long ldb = n;
			long info = 0;

			zgesv_(&n, &nrhs, dat, &lda, ipiv, gb->dat, &ldb, &info);

			m_isFactorized = true;

		};

		void SolveSVD(IComplexMatrix ^ smat, IComplexMatrix ^ umat, IComplexMatrix ^ vtmat) {

			char jobu = 'A';
			char jobvt = 'A';
			long m = m_rownum;
			long n = m_colnum;
			long lwork = Math::Max(3 * Math::Min(m, n) + Math::Max(m, n), 5 * Math::Min(m, n));
			double * rwork = new double[Math::Max(1, 5 * Math::Min(m, n))];
			doublecomplex * work = new doublecomplex[lwork];
			long lda = Math::Max(1, m);
			long ldu = m;
			long ldvt = n;
			long info;

			double * s = new double[Math::Max(1, Math::Min(m_rownum, m_colnum))];
			ComplexGeneralMatrix ^ u = (ComplexGeneralMatrix^) umat;
			ComplexGeneralMatrix ^ vt = (ComplexGeneralMatrix^) vtmat;

			zgesvd_(&jobu, &jobvt, &m, &n, dat, &lda, s, u->dat, &ldu, vt->dat, &ldvt, work, &lwork, rwork, &info);

			for (int i = 0; i < Math::Min(m_rownum, m_colnum); i++) {
				smat[i, i] = Complex(s[i], 0.0);
			}
			delete[] work;
			delete[] s;
			delete[] rwork;

			if (info == 0) {
				return;
			} else if (info < 0) {
				throw gcnew LatoolException(String::Format("the {0}-th parameter had an illegal value.", -info));
			} else {
				throw gcnew LatoolException("bidiagonal form B did not converge to zero.");
			}
		};
	};
}